Current density estimation for MEG analysis with spatial filter 
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1 Introduction 

Magnetoencephalography (MEG) [1] is applied to 
examine the higher brain functions, and there of¬ 
ten appear multiple current sources. With little 
prior information on the current source distribu¬ 
tion, the method to estimate the current sources 
is required. 

Recently, current distribution model has been 
used, and spatial filter has been expected as 
the methods. The spatial filter [2] is a signal¬ 
processing technique differing from estimation by 
inverse solution. This method produces filter 
that extracts only component of a specified source 
from measurements and inhibits the other compo¬ 
nents from the other sources. By producing the 
filter over the source region and by estimating the 
strength at a specified source, current distribution 
is obtained. 

One method to calculate the spatial filter is to use 
the temporal information from MEG data. How¬ 
ever, the time resolution that is very important 
advantage of MEG has been lost there. We sug¬ 
gested using spatial filter without temporal infor¬ 
mation. But this approach has weakness that the 
solutions are biased superficially due to the char¬ 
acteristic of the lead field of MEG sensors and 
ill-posedness of this problem. And hence, we at¬ 
tempted to improve the spatial filter. 

In this study, we proposed four methods. In first 
three methods the spatial filter itself was modi¬ 
fied, and in the last one source region was modi¬ 
fied. Each of all the methods was applied to MEG 
data, and the results were compared. Though the 
above methods did not show enough performance, 
the limitation of the spatial filter without tempo¬ 
ral information has been clarified. 

2 Basic formulation 

In this study, the physical model, a spherically 
symmetric conductor model, is adopted. Use of 


this model enables the introduction of several sim¬ 
plifications and constraints. The relationship be¬ 
tween a continuous vector current field and its in¬ 
duced magnetic field at ith MEG sensor in space 
is given by Eq.(l). 

mi = [ Gi(r') ■ J (r')dv' (1) 

Jn 

Where G is called the lead field. It describes the 
sensitivity distribution of the ith sensor and is ob¬ 
tained with the law of Biot-Savart. G is a func¬ 
tion of the relative position between the sensor 
and the source, and is inversely proportional to 
square of the distance. It depends on the con¬ 
ductivity and on the coil configuration of the sen¬ 
sor. kl is the source region. The current density 
vector is replaced by a finite number of current 
dipoles, Q k located in a three-dimensional grid 
where a dipole’s orientation and magnitude are 
determined by integrating the current field over 
the volume cell (voxel) surrounding the dipole. 
Eq.(l) becomes 

K 
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where K is the total number of dipole cells. This 
equation can be rewritten in vector-matrix form 
as the linear system. 

m = GQ (3) 

3 Spatial filter 

The spatial filter employed in this study is predi¬ 
cated upon acquisition of simultaneous MEG sig¬ 
nals from an array of sensors. The filter itself is a 
linear projection operation in which the MEG sig¬ 
nals are transformed to the source signal of a spec¬ 
ified location and a specified orientation. That is, 
it is selective for orientation of dipole as well as 
for position. 



A specified location and orientation of the &th 
source is regarded, and it’s filter is w * The filter 
output for &th source is: 

M 

s k = J2 W ki m i = w k m ( 4 ) 

Eq.(4) of all sources are calculated, and the cur¬ 
rent density distribution is obtained over the 
source region. 

These filters are obtained as bellow. Prom Eq.(3) 
and Eq.(4), 

s k = wlGQ (5) 

The properties of this filter are that it extracts 
only components of &th source from measure¬ 
ments, and that it inhibits components arriving 
from the other sources. Thus, 

W k Gj fijk 

where 8 k is defined as having a value of one for k 
and a value of zero everywhere else. And |<S| = 1. 
The filter coefficients are computed to minimize 
the difference in a least-squares sense. 

E p r w k 8 k \' 2 (6) 

The coefficients minimizing F are found at 
dF/dw k = 0. The solution for a useful array 
of the spatial filter coefficients may be found: 

w k = {GG T )~ 1 G8 k = ( GG T )~ l b k (7) 

where b k = G8 k is the forward solution for the 
physical sensors due to a current dipole, 8 k , of 
unit strength at the &th source position. 


4 Improvements of spatial filter 

The solutions by the original spatial filter are bi¬ 
ased superficially due to the characteristic of the 
lead field of magnetic sensors. Thus we propose 
the methods to modify the spatial filter as follows. 

4.1 Method 1 

The spatial filter tuned for the source at a speci¬ 
fied location and an orientation is derived so as to 
be sensitive. Thus, if the spatial filter is applied 
to the forward solution, b, due to a current dipole 
at the &th source position, the solution becomes 
ideally the source itself, that is, w T b —» 1. By 
adding the constraint to Eq. (6), the spatial filter 
w is modified. 


F = | \G r w k -8 k || 2 4r-J| \blw k - l|| 2 (8) 

w k = (1 + M) (GG T + X^kbiy 1 b k (9) 

Where the parameter Ai controls the influence of 
the constraint in Eq.(8). 

4.2 Method 2 

In addition to method 1, the spatial filter is modi¬ 
fied to inhibit components arriving from the other 
sources. And it is weighted with 1/r 2 to correct 
the depth of the solutions, r is the distance be¬ 
tween origin and source position. The strength of 
G multiplied by 8 k in Eq.(7) was considered to be 
the factor that the more superficial sources were 
estimated larger. So, 1/r 2 has been selected as 
the weight, because the lead field is inversely pro¬ 
portional to square of the distance between the 
source and the sensor. The spatial filter is modi¬ 
fied. 

F = ^G T w k - ±8 k || 2 | B T w k - 8' k f (10) 

w k = (1 + A 2 ) (GG T + A 2 BB t ) _1 b k (11) 

Where B is the matrix of the forward solutions 
due to a current dipole of unit strength at all po¬ 
sitions and all orientations. And the delta func¬ 
tion 8' k is defined as having a value of one for k 
and a value of zero everywhere else. 


4.3 Method 3 

The entropy for power distribution of the current 
density is used as the regularization term. By 
minimizing the entropy, pj, the power of current 
distribution and in the range of 0 1, approaches 
closer end value gradually. Thus, the current dis¬ 
tribution becomes localized further. 

2 K 

F = | G T w k - £fc|| - 1( W-r ( 12 ) 


= ( w l b i) 
Ef«b .) 2 


(13) 


It is impossible to obtain these filter coefficients 
w by finding at dF/dw = 0 as in the above same 
ways. Thus, for minimization of Eq.(12) the con¬ 
jugate gradient method is used, and then w is 
obtained. 




4.4 Method 4 


The source region is gradually reduced, and any 
influences on current density distribution esti¬ 
mated with original method are evaluated. By 
detecting the change between the cases that the 
current source is inside and outside the evaluated 
source region, the depth of the current source is 
determined. The summation of the squared resid¬ 
uals between measurements and estimated mag¬ 
netic fields, i.e. the error of magnetic fields, is 
used as an indicator for detecting this change. 


e = 100 


HiLiimj -rhj ) 2 
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(14) 


Where m* is the measurement of ?'th sensor and 
rhi is the value calculated by estimated current 
distribution. 


5 Simulation 

To compare the proposed method with the orig¬ 
inal method, three simulations were performed. 
The Whole-Head Type 64 channels MEG System 
was used in this study. It was assumed that the 
source region was hemisphere with the radius of 
7 cm and were digitized with the intervals of 1 
centimeter between sources. Each methods were 
applied to MEG signals simulated with Eq.(l). 
Measured noises were added to those signals. 

5.1 Simulation 1 

Method 1 and 2 were compared with original one. 
A current source was assu m ed to be at (0,0,5). 
The SNR was 11.83 dB. The parameter Ai was 
332 in method 1 and A 2 was 5.60 x 10 -3 in method 
2 . 

The results are shown in Fig.l. In the original 
method the solution were became biased superfi¬ 
cially. In the method 1, the current distribution 
became maximized at the actual source position, 
(0,0,5). However, region of high strength of cur¬ 
rent distribution was expanded in all directions. 
For the constraint — l|| that the outputs 

approach ideal values the ability of the spatial fil¬ 
ter to inhibit the components arriving from other 
sources has been lost. 

In the method 2, the correction of the depth in 
Eq.(10) made partly the region of high current 
distribution deeper in the radial direction, but the 


current distribution became maximized at (0,0,7) 
as well as original one. The term of the constraint 
is almost the same as — 8^ || in the spa¬ 

tial filter and is inefficacious to improve the esti¬ 
mation of the depth. Next, to correct the depth 
of the solutions, the weight 1/r 2 was used in the 
consideration of the lead field. However, the ill- 
posedness of this problem would greatly affect on 
the weakness and the weight did not show enough 
performance. 



a) original method b) method 1 



c) method 2 


Figure 1: Current density distributions estimated 
with each of the methods. They are displayed as 
arrow maps of cross sections of current distribu¬ 
tions in the xz plane with y=0. Each black arrow 
denotes the estimated current density at the posi¬ 
tion and a gray arrow denotes the actual current 
source. 


5.2 Simulation 2 

Method 3 was compared with original one. How¬ 
ever, since the spatial filter regularized by entropy 
took very long time to be computed, it was as¬ 
sumed that the interval between sources was 2 
cm. The number of sources became 107. A cur¬ 
rent source was assumed to be at (0,0,4). The 
SNR was 8.85 dB. 

The results are shown in Fig.2. Both the results 
were almost same. However, the method 3 was 
planed so as to improve the localization of the 
current distribution. So, the estimated source 
strengths were compared. The results are shown 
in Fig.3. The difference was minus near maximum 
source. Thus it is said that the current distribu¬ 
tion have been more localized by the regulariza¬ 
tion term by entropy. 
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original method method 3 

Figure 2: Current density distributions estimated 
with each of the methods. They are displayed as 
arrow maps of cross sections of current distribu¬ 
tions in the x-z plane with y=0. 



Figure 3: Plot of estimated source strengths in 
ordinate and sources in abscissa. Sources are ar¬ 
ranged by distance from the actual source starting 
with the nearest one. Estimated source strengths 
are value that the source strengths with original 
method are subtracted from ones with method 3. 

5.3 Simulation 3 

Simulation was performed to investigate method 
4. The original method is used there. The radius 
of the source region was changed from 7.0 to 4.0 
cm. However, the intervals between the sources 
were assumed to 0.5 cm to evaluate exactly the in¬ 
fluence. A current source was assumed to (0,0,5.8) 
and the SNR was 14.60 dB. 

The result of one current source is shown in Fig.4. 
From Fig.4, the error changed in the vicinity of 
the dotted line. Thus, it is said that the depth 
of the current source was estimated by detect¬ 
ing the change of Eq.14. The current distribu¬ 
tions of the radius 6.0 and 5.5 cm are shown in 
Fig.5. Thus, the position could be estimated to 
be around (0,0,6) that the current density is max¬ 
imized in Fig.5. However, method 4 could be not 
enough to estimate in the existence of some cur¬ 
rent sources. 


Figure 4: Error of magnetic fields in Eq.lf with 
radius of the source region. The dotted line indi¬ 
cates the position of the actual current source. 



a) 6.0 cm b) 5.5 cm 


Figure 5: Current density distributions projected 
onto xz plane. The only arrows over the half value 
of maximum current distribution are shown. 


6 Conclusion 

In this study, we attempted to improve the esti¬ 
mation of the depth of current sources by modi¬ 
fying the spatial filter without temporal informa¬ 
tion. As a result, though the above methods did 
not show enough performance, the limitation of 
this approach has been clarified. 
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